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We present a fast and accurate method to calculate vibrational properties for mechanically un- 
stable high temperature phases that suffer from imaginary frequencies at zero temperature. The 
method is based on standard finite-difference calculations with optimized large displacements and 
is significantly more efficient than other methods. We demonstrate its application for calculation 
of phonon dispersion relations, free energies, phase transition temperatures, and vacancy formation 
energies for body-centered cubic high-temperature phases of Ti, Zr, and Hf. 



I. Introduction 

Hexagonal close packed (hep) to body-centered cubic 
(bee) phase transitions are characteristic of a number 
of transition and rare earth metals, including titanium, 
zirconium, and hafnium. Although room-temperature 
applications of these materials are limited by the ex- 
pense of processing them, their high-temperature prop- 
erties make them ideal candidates as the basis of high- 
temperature alloys. However, only limited experi- 
mental datai~— have been collected on the bee high- 
temperature phase of these systems, and they are inac- 
cessible to traditional zero-temperature ab-initio calcula- 
tions due to their mechanical instability— Although the 
recently proposed self-consistent ab-initio lattice dynam- 
ics (SCAILD)^ and molecular dynamics (MD) based^ ap- 
proaches have proved successful in calculating phonon 
dispersions for these materials, convergence issues with 
free energy and phonon slopes in the long-wavelength 
limits prevent accurate calculation^ of other phase prop- 
erties relevant to high-temperature applications (e.g. 
elastic constants and formation energies of vacancies). 
As a result, the computational cost of these advanced 
calculation procedures sets limits to their applicability, 
and a supplemental calculation method is necessary for 
complete phase analysis. Thus, there is significant mo- 
tivation to develop a simplified method of dealing with 
mechanically unstable systems of potentially large size 
that quickly produces reasonable ab-initio results. 

Current ab-initio methods to deal with mechanical in- 
stability in these systems have focused on a simulated 
MD approach by displacing atoms around their equilib- 
rium positions and using the resulting structures to cal- 
culate dynamical matrices for phonon calculations. The 
SCAILD method proposed by Souvatzis et alii is lim- 
ited by its high symmetry requirement for reasonable 
calculation convergence, and is thus inefficient at pro- 
ducing calculations beyond phonon dispersions for per- 
fect structures^ A second calculation method proposed 
by Hcllman et al. utilizes an MD approach to extract dy- 
namical matrices, but is computationally intensive due 
to the large number of MD steps it necessitates for 
convergence^ Due to their computational demands, nei- 
ther method has been used to date to calculate the ener- 
getics of defects in mechanically unstable systems. 

In this letter, we propose a method to calculate ac- 



curate thermodynamic properties of mechanically un- 
stable high-temperature phases using ab-initio theory. 
This method utilizes large-displacement phonon calcula- 
tions with a displacement determined by analysis of the 
phonon amplitude-energy curve. Because these calcula- 
tions arc computationally very economical, more com- 
plex calculations of defect energies, alloying effects and 
transition temperatures are possible than using previous 
ab-initio methods. Furthermore, we find for the high- 
temperature bec phases that we have studied that the 
agreement of phonon dispersions and phase transition 
temperatures with experimental values is at least as good 
as with previously published ab-initio methods. 

In high-temperature bee phases, where the mechan- 
ical instability occurs as a result of decreased energy 
when a (£,—£,0,0,0,0) strain is applied to the calcula- 
tion cell^, the energy-strain response in this direction 
can be approximated by a quartic function with negative 
curvature at zero strain. The corresponding elastic con- 
stant c', given by that curvature, is thus negative. Due 
to the relationship between the slopes of the (acoustic) 
phonon branches and the elastic constants, especially for 
the transverse branch in [110] direction which is given by 
(Cn — C12) /2, the negative elastic constant corresponds 
to imaginary frequencies in the phonon dispersion, which 
is determined from the basic assumption that the energy 
vs. atomic displacement curves have positive curvature 
(minimum at the equilibrium position) and are parabolic 
(harmonic). Phonon frequencies and eigenvectors are de- 
termined by the curvature of the energy vs. atomic dis- 
placement curve at the equilibrium position and can be 
determined by density-functional perturbation theory^ or 
small finite displacements.— 

The appearance of imaginary frequencies does not only 
interfere with determining phonons per se, but also with 
free-energy based calculations, which for solids are eas- 
iest obtained within the quasiharmonic approximation. 
There, the free energy contribution to the total en- 
ergy is calculated from the phonon density of states.— A 
small-displacement phonon calculation of the free energy 
contribution to mechanically unstable structures is thus 
flawed through its inclusion of nonphysical frequencies. 

The fundamental question at this point is why the 
high-temperature phases, despite being mechanically un- 
stable, become stable at sufficiently high temperatures. 



As has been discussed previously, this is due to the fact 
that the vibrational amplitudes of the atoms increase 
with temperature. This has two consequences: firstly, 
an atom no longer moves in the potential landscape of 
the other atoms at or close to their equilibrium lattice 
positions, but displaced from them (phonon-phonon cou- 
pling). Secondly, the atom explores a larger part of the 
energy landscape through its higher vibrational ampli- 
tude (phonon anharmonicity). If the second effect dom- 
inated the first, a calculation that includes these anhar- 
monic effects could be based on an improved displace- 
ment strategy around the equilibrium positions of the 
atoms. This would avoid the major computational over- 
head in the methods proposed to date^, which consist 
of relating displacements for the other atoms in the sys- 
tem to the atom displaced for a phonon calculation. In- 
deed, such a strategy may be successful and not ham- 
pered by the instability in direct vicinity of the equilib- 
rium position considering that an oscillating atom has its 
largest velocity when it passes through the equilibrium 
site, while its velocity decreases to zero at the turnaround 
points, where it thus spends the majority of its time. 
Therefore, we propose to use atomic displacements for 
the phonon calculations large enough to probe the cur- 
vature of the energy landscape far away from the equi- 
librium position. 

Indeed, at large displacements, the quartic term of 
the energy-strain function eclipses the negative quadratic 
term, leading to an approximately constant positive cur- 
vature at sufficient distance (e.g. at ^0.5 A in Fig. [TJ). 
Since a constant curvature corresponds to harmonic mo- 
tion and enables the phonon concept, the shape of the 
curve can thus be approximated by only using a quadratic 
term with little loss of detail. This suggests that a 
phonon calculation of the free energy contribution us- 
ing a large displacement in the quartic-dominated regime 
could be used with the quasi-harmonic approximation to 
provide the free energy contribution to the total struc- 
ture energy, and that the properties of a structure thus 
stabilized could be calculated using ab-initio methods. 
By comparing the results of a large-displacement calcu- 
lation to experiment, one then can also probe in a simple 
way if the dominating anharmonic effect comes from the 
large atomic displacements (with lesser effect from what 
the surrounding atoms do), or vice versa. 

II. Methodology 

To calculate the approximate large displacement nec- 
essary to force the phonon calculation into the quartic- 
dominated regime, VASP— was used on a frozen-phonon 
type supercell for A-point phonons 1 ^ to calculate the in- 
ternal energy of the structure as a function of phonon 
displacement^ The A-point phonon was chosen since 
the transverse mode at the A-point is the primary marker 
of mechanical instability in bec transition metalsAi For 
these calculations, we used projector-augmented wave 
potential o 13 ' 14 with a 10 x 10 x 10 Monkhorst-Pack k-point 
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FIG. 1. The (a) energy- vs. -phonon amplitude curve and its 
(b) second derivative for an iV-point frozen phonon cell^ dis- 
placement of the bec cell of Hf (solid) and Mo (dashed). 



mesh. Calculations were performed for bec titanium, zir- 
conium and hafnium, as well as several mechanically sta- 
ble bec metals. 

Figure [JJa) shows the calculated internal energy as a 
function of the A-point phonon amplitude in the cho- 
sen displacement direction for the Hf and Mo systems. 
To determine the appropriate phonon displacement for 
large-displacement calculations, second numerical deriva- 
tives of this internal energy data are calculated using 
a central-difference method, shown in Figure 1(b). In 
the case of Hf, these calculations give an indication of 
the point at which the amplitude-energy curve enters the 
quasi-parabolic high-temperature regime (in this case be- 
tween ~0.4 and 0.7 A, where the second derivative be- 
comes approximately constant), allowing harmonic- type 
phonon calculations and use of the quasi-harmonic ap- 
proximation to determine the system's free energy. By 
contrast, the second derivative for Mo is maximally posi- 
tive at values around zero displacement, which therefore 
defines its quasiharmonic regime. For automatic detec- 
tion of an appropriate displacement value, we found that 
the maximum value of the second derivative identifies 
well the necessary displacement to calculate a realistic 
phonon dispersion. This interpretation matches well with 
phonon approaches for stable materials, which exhibit 
maximum second derivative values at minimal displace- 
ments. Alternatively, one can select the displacements 
just large enough that all the imaginary frequencies from 
the phonon dispersions are eliminated. Typically, the last 
frequencies to turn real are located not at the A-point, 
but close to the T-point. 

Once the necessary displacement magnitude was deter- 
mined for each system by eliminating all the imaginary 



TABLE I. Calculated vs. experimental values (in paren- 
theses) 043, I15HI7I ] of phase transition temperature T tra ns, 
bulk modulus K, and lattice parameter a for the bcc high- 
temperature phases of Ti, Zr, and Hf. Also given are large- 
displacement values determined from maximum values of the 
second derivative of TV-point frozen- phonon cells (u2d) and 
from a no- imaginary- frequency condition (u n i) used for these 
calculations. 





U2d (A) 


Unl (A) 


Ttrans (K) a (A) B (GPa) 


Ti 


0.60 


0.88 


1200 (1155) 3.276 (3.33) 87.1 (87.7) 


Zr 


0.57 


0.54 


1050 (1136) 3.557 (3.551) 94.8 (96.7) 


Hf 


0.63 


0.68 


2100 (2016) 3.521 (3.625) 119 (112.3) 



frequencies from the phonon dispersion, force constants 
and dynamical matrices were determined from finite- 
displacement VASP calculations within 4x4x4 supercells 
of the bcc primitive cells using the PHON code.— The cal- 
culations were performed with 6x6x6 Monkhorst-Pack 
k-point meshes for a range of lattice parameters. Us- 
ing the quasi-harmonic approximation in PHON^, these 
calculations were used to determine the vibrational con- 
tribution to the system's free energy over a temperature 
range that included the experimental transition tempera- 
ture. For all materials considered, an additional calcula- 
tion of the phonon density of states was conducted on re- 
laxed supercells composed of 3 x 3 x 2 hexagonal unit cells 
with a r-centcred 6x6x9 Monkhorst-Pack k-point mesh 
using a 0.03 A phonon amplitude. The lattice parame- 
ter was then determined by minimizing the free energy 
with respect to volume. The minimum of the free energy 
vs. volume curve at each temperature was then used to 
determine the phase transition temperature. The bulk 
modulus was calculated using a Birch-Murnaghan fit for 
temperatures across the literature range of phase stabil- 
ity and plotted to determine its temperature dependence. 

To calculate the free energy of formation of a vacancy 
in the metallic bcc system, F f (T) = F V {T) - ^-F p {T), 
the free energy as a function of temperature for 3x3x3 
54-atom bcc super cells with (V) and without (p) va- 
cancies was calculated using a 6x6x6 Monkhorst-Pack 
k-point mesh and the large-displacement values from the 
perfect cells. For this demonstration here, the vacancy 
structures were not relaxed. 

Table |T] shows the calculated transition temperature, 
lattice parameter, and bulk modulus obtained using the 
large displacement phonon approach wherein no imagi- 
nary frequencies were present in the phonon dispersion, 
as well as experimental values. In addition, the displace- 
ments are displayed at which no imaginary frequencies 
appear in the phonon dispersion, and at which the sec- 
ond derivative of the phonon amplitude-energy curve was 
maximal. 

III. Results 

For bcc titanium, the first displacement at which no 
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FIG. 2. Large-displacement calculated phonon dispersions for 
(a) Ti, (b) Zr, and (c) Hf vs. experimental phonon energies 
from [J^. Solid lines are calculated with 4x4x4 phonon 
supercells, the dashed line for Hf with 8x8x8. 



imaginary frequencies were observed in the phonon dis- 
persion was 0.88 A; this value did not correspond to the 
absolute maximum of the second derivative, which oc- 
curred at approximately 0.6 A, but did correspond to 
a second local maximum of the second derivative which 
occurred at approximately 0.85 A. The transition tem- 
perature was calculated to be 1200 K, which is within 
4% of the experimental value of 1155 K— , and closer 
to experiment than that calculated using the SCAILD 
method, 1107 K— . The lattice parameter was calculated 
to be 3.234 A at K and 3.276 A at the calculated transi- 
tion temperature, while the experimental lattice param- 
eter at this temperature was 3.33 A. The coefficient of 
thermal expansion was calculated to be 6.7 x 10~ 6 K -1 , 
compared to an experimental coefficient of thermal ex- 
pansion for the BCC phase of 8.6 x 10~ 6 K" 1 ^ The 
calculated phonon dispersion using the 0.88 A phonon 
amplitude is displayed in Figure [2ja), with experimental 
data points from Ref. [l|. 

For bcc zirconium, the first displacement at which no 
imaginary frequencies were observed in the phonon dis- 
persion was 0.54 A; this value agreed well with the cal- 
culated maximum of the second derivative, which was 
approximately 0.57 A. The transition temperature was 
calculated to be 1050 K, which is 8% lower than the 
experimental value of 1136 K— , whereas the transition 
temperature calculated using the SCAILD method, 916 
K-£, is more than 19% too low. The lattice parameter 
was calculated to be 3.557 A at the calculated transition 
temperature and 3.562 A at 1150 K, while the exper- 
imental lattice parameter near the transition tempera- 
ture was 3.551 A. The coefficient of thermal expansion 
was calculated to be 1.34 x 10~ 5 K _1 (no corresponding- 
experimental measurements were found). The calculated 
phonon dispersion using the 0.54 A phonon amplitude is 



displayed in Figure [2](b) . in comparison to experimental 
data points.— 

For bec hafnium, the first displacement at which no 
imaginary frequencies were observed in the phonon dis- 
persion was 0.68 A; this value agreed well with the calcu- 
lated maximum of the second derivative, which was ap- 
proximately 0.63 A. The transition temperature was cal- 
culated to be 2100 K, which is 4% higher than the exper- 
imental value of 2016 K— , while the transition temper- 
ature calculated using the SCAILD method, 1660 K— , is 
18% too low. The lattice parameter at 2400 K was calcu- 
lated to be 3.521 A, which is 3% smaller than the experi- 
mental value of 3.625 A at that temperature. The coeffi- 
cient of thermal expansion was calculated to be 7.9 x 10 -6 
K _1 , compared to an experimental coefficient of thermal 
expansion for the BCC phase of 11 x 10" 6 K -1 ^ The 
calculated phonon dispersion using the 0.68 A phonon 
amplitude is displayed in Figure EJJc), with experimental 
data points^. Additional phonon calculations were per- 
formed on the Hf bec system using larger supcrcells to de- 
termine whether increases in cell size would improve the 
agreement of results with experimental values. Within 
the feasible limits of computational resources available 
for this study, no significantly better agreement was ob- 
served from increasing the calculation cell size past a 
supercell composed of 4x4x4 bec primitive cells. The 
results of expanding the supercell size from 4x4x4 to 
8x8x8 bec primitive cells is shown in Fig. [2](c) . 

The temperature dependence of bulk moduli for Ti, 
Zr, and Hf was calculated using the large-displacement 
method with Birch-Murnaghan fitting across a tempera- 
ture range as previously described. The three tempera- 
ture dependencies, displayed as normalized to each ma- 
terial's transition temperature in Figure [3j exhibit nearly 
perfect quadratic relationships. Also plotted in Figure [3] 
are experimental data for the bulk moduli of the three 
materials, shown as solid points. Calculated values for 
Ti and Zr exhibit the typically expected decrease in bulk 
modulus with increased temperature. The experimen- 
tal data for Ti (Fig. [3ja)) show a large scattering and 
indicate significant error bars for the measured values, 
which are not reported in the experimental source fl9| . 
Within these assumed error bars, our temperature de- 
pendence is in agreement for Ti and Zr, where values 
have been measured for more than one temperature. Cal- 
culated values for Hf indicate an anomalous increase in 
bulk modulus with increased temperature and show good 
agreement with the one experimental value available Q 
for a temperature of 2073K. While it is possible that this 
represents a calculation artifact or results from the Birch- 
Murnaghan fits used, if physically valid it would indicate 
an internal structural change that may be related to the 
mechanism of the HCP-BCC phase transition. 

The free energy of formation for unrelaxed vacancies 
in bee Ti, Zr, and Hf were calculated using the large- 
displacement method with the displacements previously 
described. By fitting F f (T) = E f - TS>1£ in the bec 
stability range, the energy Ef and entropy Sf of forma- 
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FIG. 3. Calculated temperature dependence of bulk modulus 
for (a) Ti, (b) Zr, and (c) Hf, normalized to the transition 
temperature of each metal, with corresponding experimental 
values from @,H and El. 



tion of the vacancy were determined. In Ti, the vacancy 
energy of formation Ef was 2.05 eV, and the entropy 
of formation Sf was 8.15/cb. In Zr, the vacancy energy 
of formation Ef was 2.06 eV, and the entropy of forma- 
tion Sf was 2.76/cb. In Hf, we found Ef = 2.20 eV and 
Sf = 0.89/cs. The values for the formation energy deter- 
mined in this way are smaller by 0.9%, 0.6%, and 5.9% 
than the zero-temperature formation energies of 2.07, 
2.08, and 2.34 eV for Ti, Zr, and Hf respectively. While 
the agreement between fit and zero-temperature forma- 
tion energy for Ti and Zr is similar to more harmonic bee 
metals, where we calculate a decrease of 0.2%, 0.4%, and 
1.0% for the examples of W, Mo, and Cr, respectively^! 
Hf shows a strong discrepancy coming from the fact that 
its phonon density of states shows a strong decrease in 
the high-frequency regime upon introduction of a vacancy 
which is not observed for Ti and Zr. Additionally, it 
would not be possible to calculate the entropy of vacancy 
formation without a reliable temperature dependence of 
free energy derived from a converged phonon density of 
states, as calculated using this method. We find that the 
entropy of formation scales inversely with the melting 
temperature of the three metals, which is sensible, since 
the melting temperature is a measure for the influence 
of the lattice vibrations on the stability of the lattice the 
higher the melting temperature, the smaller the effect 
of lattice vibrations on the order of the atoms with re- 
spect to their equilibrium position, and thus the smaller 
the influence of vibrational disorder on the formation of 
vacancies. In order to reliably calculate the relative mag- 
nitude of these contributions for anharmonic structures 
at high temperatures, full free-energy calculations are a 
necessity, and our results show that they are well within 
the capabilities of the method proposed here. 

IV. Conclusions 
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In this paper, the large-displacement phonon ap- 
proach has been shown to produce results in good agree- 
ment with experimental values for a number of high- 
temperature metallic systems. However, it cannot be 
expected to be applicable to all systems in the same 
way presented here, as the causes underlying mechani- 
cal instability in these systems may differ, although pre- 
liminary calculations for other systems up to now look 
promising. Our methodology can only be applied to those 
systems in which a parabolic regime can be reached by in- 
creased phonon amplitude, as determined in the displace- 
ment calculation steps. For these systems, it produces 
results that have excellent agreement with experimental 
values at a relatively low computational cost. For the 
other systems, large-displacement calculations can help 
to easily show that phonon-phonon interactions domi- 
nate the effect of a locally anharmonic potential. 

In summary, the primary advantages in the large- 
displacement method for calculating the energetics of me- 
chanically unstable high temperature phases and defects 
are its speed and simplicity. Due to the computational 
demands of both the SCAILD£ and MD& approaches, 



the cell volume is at best estimated from an extrapola- 
tion from a K DFT relaxation, which can lead to large 
errors as shown here. In contrast, the large-displacement 
method provides a fast self-consistent approach that al- 
lows phonon calculations on a number of different cal- 
culation cell volumes simultaneously, allowing for true 
incorporation of temperature effects due to lattice ex- 
pansion and atomic vibration to yield more accurate val- 
ues for thermodynamically determined quantities such 
as phase transition temperatures. The incorporation of 
these effects via the large-displacement method seems to 
allow for results in the same range of accuracy as those 
obtained for mechanically stable phases. Additionally, 
the significantly reduced computational cost allows for 
more complex calculations to be performed, including de- 
fect calculations, that were previously beyond the reach 
of ab-initio methods. 
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